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I.  Introduction 


The  Air  Force’s  Engine  Structural  Integrity  Program  (ENSIP)  (Ref  1) 
requires  determination  of  damage  tolerance  for  jet  engine  components  in  order 
to  allow  more  economical  rejection  criteria  to  be  adopted.  To  this  end,  means 
have  been  developed  for  predicting  fatigue  crack  growth  in  jet  engine  com¬ 
ponents  such  as  turbine  disks  made  of  nickel-based  superalloys  and  operating  at 
elevated  temperatures.  The  presence  of  time-dependent  plastic  deformation 
greatly  affects  crack  propagation  rates,  particularly  at  elevated  temperatures 
and  thus  must  be  account,  d  for  when  modelling  crack  growth  in  turbine  materi¬ 
als.  - - 


This  plastic  deformation  in  such  materials  takes  two  forms:  That  due  to 
viscoplasticity  and  that  due  to  creep.  Viscoplasticity,  as  defined  herein,  occurs 
at  high  stresses  and  is  characterized  by  large  plastic  deformation  rates  accom¬ 
panied  by  changes  in  the  material’s  internal  state,  typified  by  strain  hardening 
or  softening  in  which  the  material  becomes  more  or  less  resistant  to  plastic  flow. 
Creep,  as  defined  herein,  occurs  at  lower  stresses  and  is  characterized  by  small 
plastic  strain  rates  while  the  material’s  internal  state  remains  essentially 
unchanged.  Viscoplastic  deformation  rate  and  creep  rate  are  both  functions  of 
the  magnitude  and  frequency  of  the  applied  strains  and  stresses.  Finite  element 
methods  (FEM)  can  be  used  to  predict  crack  growth  in  a  component,  but  must 
incorporate  a  constitutive  model  that  accounts  for  the  viscoplasticity  and  creep 
which  will  occur  at  notches,  crack  tips,  and  other  stress  concentrations  when 
the  applied  stress  is  sufficiently  high. 

A  number  of  methods  for  modelling  viscoplasticity  and  creep  exist,  some  of 
which  calculate  vise  op  las  t  iv  it  v  and  creep  via  separate  equations  while  others 
unify  viscoplasticity  and  c  i  >■  <  p  with  a  single  equation  which  is  valid  fe.i  high  and 

In  i'  ality  the  transition  between  creep  and  viseoplas- 


low  inelastic  strain  rates 


ticity  often  occurs  gradually  over  a  range  of  stresses  and  strains,  thereby  mak¬ 
ing  a  unified  model  more  consistent  with  actual  material  behavior.  For  this  rea¬ 
son,  a  unified  model  was  sought.  One  such  model  which  the  Air  Force  has  inves¬ 
tigated  for  crack  growth  prediction  is  the  Bodner-Partom  constitutive  law  (Ref 


2)- 


— ^  The  purpose  of  this  thesis  was  to  investigate  the  frequency  response 
aspects  of  the  Bodner-Partom  constitutive  law's  behavior  and  to  compare  its 
results  with  those  of  other  models  and  to  cyclic  and  non-cyclic  uniaxial  tensile 
test  data. 


The  investigation  proceeded  in  several  stages.  First,  two  FEM  codes  were 
developed  to  model  uniaxial  viscoplasticity  and  creep,  one  using  Bodner-Partom 
constitutive  theory,  the  other  a  combination  of  Overstress  Law  and  Norton’s  Law 
for  Secondary  Creep.  The  codes  were  then  validated  by  comparison  of  their 
predictions  with  those  of  other,  proven  codes  and  with  data  taken  from  uniaxial 
tensile  tests  upon  specimens  subjected  to  cyclic  stress  and  strain  over  a  range 
of  frequencies  and  amplitudes.  The  programs  were  then  used  to  extrapolate 
uniaxial  material  behavior  for  conditions  not  covered  by  test  data  and  to  show- 
material  behavior  at  a  stress  concentration  while  under  cyclic  stress  using  a 
simple  three-bar  structural  model. 
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II.  Literature  Search 


Available  literature  was  reviewed  in  order  to  identify  other  ways  by  which 
viscoplasticity  has  been  modelled  and  to  find  which  models  have  been  most  use¬ 
ful  for  crack  growth  prediction. 

The  review  turned  up  a  number  of  viscoplasticity  models,  all  of  which  have 
two  elements  in  common:  a  means  of  relating  inelastic  strain  rate  to  applied 
stress,  and  a  means  of  accounting  for  strain  hardening  as  a  function  of  accumu¬ 
lated  inelastic  strain.  Typical  viscoplastic  strain  rate  equations  include: 


Power  Law  (Refs  3,4,5): 


ip  =  C 


a 

Y 


Overstress  (Refs  5,6): 

tr-X 


=  C 


ip  =0 


(  *  >X) 

(  <7  <X) 


Exponential  (Refs  2, 4, 7, 8): 

(f)- 


ep  =  C  e 


Where  a  is  the  applied  stress,  C  and  n  are  material  constants,  Xdetermines 
the  onset  of  plasticity  for  the  Overstress  model,  and  Z  is  a  state  parameter  used 
to  show  hardening.  Isotropic  or  non-direc  tional  hardening  can  be  account  ed  for 
by  making  Z  a  function  of  accumulated  plastic  strain  (Refs  2, 3, 5, 7),  while 
kinematic  or  directional  hardening  can  be  portrayed  by  adding  to  Z,  terms 
whose  value  depends  upon  the  direction  of  the  applied  stress  (Refs  5,7,8).  Titer- 


III.  Description  of  Methods 


Introduction 

In  this  section,  it  will  be  shown  how  inelastic  strain  rates  are  calculated  and 
how  hardening  and  recovery  are  portrayed  by  the  Bodner-Partom  model  and  by 
the  combined  Ove  rstress/ Norton ’s  Law  model.  In  addition,  time  step  selection 
and  other  numerical  stability  criteria  will  be  described. 

The  Bodne  r-Partom  Constitutive  Law 

Bodner  and  Partom  developed  a  constitutive  relationship  between  applied 
stress  and  the  rate  of  plastic  deformation  (Ref  2).  The  Bodner-Partom  model 
can  be  used  to  predict  viscoplasticity  and  creep  and  includes  the  effects  of 
strain  hardening/softening  and  thermal  recovery. 

The  Bodner-Partom  model  relates  plastic  deformation  rate  to  applied  stress 
by  using  a  state  variable  Z  in  the  following  manner: 

*  ipiiepij  =  D2p(J2,Z)  (1) 

where : 

-  D 2P  is  the  second  invariant  of  the  plastic  strain  rate  and  is 
give n  by 

D/.Cl.x  p(-(^i-)'"±L  (2) 

-  J 2  is  the  second  invariant  of  the  deviatoric  stress. 

-  D0  is  the  limiting  value  for  the  strain  rate  (for  small  strain 


rates,  generally  set  at  an  arbitrarily  high  number  such  as 
I.0X104  ser) 


*  the  parameter  n  determines  the  yield  stress  at  a  given  strain 
rate . 


(3) 

(4) 


Isotropic  strain  hardening  and  softening  is  modelled  by  allowing  the  param¬ 
eter  Z  to  evolve  as  plastic  strain  accumulates: 


Z=Z1-(Zl-Z0)ex  p(-m  Wp  ) 


(5) 


where  Wp  is  the  plastic  work  and  m  is  a  material  constant  that  determines  the 
rate  of  change  of  Z  and,  by  extension,  the  shape  of  the  stress-strain  curve  for  a 
given  stress  or  strain  rate.  Z0  represents  the  initial  value  of  Z,  and  Z j 
represents  the  limiting  value  which  Z  will  approach  as  plastic  strain  occurs  .If  the 
material  exhibits  strain  hardening,  Zx  is  greater  than  Z0.  If  the  material  exhi¬ 
bits  strain  softening,  Zx  is  less  than  Z0. 

In  tnis  investigation,  the  plastic  work  %  is  defined  to  include  the  net  plastic 
strain  energy  plus  an  additional  term  to  account  for  thermal  recovery  of  har¬ 
dening  at  elevated  temperatures: 


0  0 


(uniaxial  case' 


where 


where  Z,  represents  a  lower  limit  of  Z,  and  A  and  r  are  additional  material  con¬ 


stants  . 

Viscoplastic  stress  and  strain  are  calculated  by  Euler  numerical  integration 
over  time.  During  each  time  step  i,  Equations  (3)  through  (8)  are  performed  for 
each  element: 


Z*-Z,-(ZI-Z0)exp[-m  Wj 

[D/]’=Do2exP  L(  !—)*»(  “±i-) 


n 


[  < p  **  ] '  = 

[dt?x]i={ip  „y  dr 


'vtr'+*iAd€?,\ ,+(z,»]< 


m^-Z1) 


(9) 

(10) 

(11) 

(12) 

(13) 


The  material  constants  are  determined  via  test  data  for  the  required  tem¬ 
perature.  Stouffer  (Ref  IV)  shows  how  to  obtain  material  parameters  Zx,  Z0l  Z, , 
A,  m,  n,  and  r  from  uniaxial  tensile  stress-strain  tests  and  creep  tests  at 
different  stress  and  strain  levels.  The  specific  material  parameter  values  used 
in  this  investigation  are  shown  in  Table  I  for  turbine  materials  IN  100  and  Inconel 
718. 

Note  that  for  this  investigation  only  isotropic  or  non-direc tional  hardening 
has  been  incorporated  into  the  Bodne r-Par tom  model.  Kinematic  or  directional 
hardening  requires  the  calculation  of  several  additional  material  parameters 
and  their  evolution  rates  during  each  time  step  and  so,  for  simplicity,  was  not 
modelled.  If  required,  the  methods  described  by  Lindholm,  et  al  (Ref  7)  and  by 
Beaman  (Ref  8)  could  be  incorporated  into  the  model 


Overstress  Theory  and  Norton's  Law  for  Secondary  Creep 

The  overstress  formulation  used  was  that  developed  by  Perzyna  (Ref  12).  In 
this  case,  viscoplastic  strain  rate  is  determined  as  a  function  of  applied  stress 
by  assuming  the  following  relationship: 


i'iy-7 


V377 


ay 


Sit 


(2~  \/3  J  2  j 


=  0 


( \/3  J  2  >  a  y) 
( V3 J 2  <<*y) 


(14) 


where  q  and  n  are  material  constants  determined  from  uniaxial  tensile  stress- 
strain  test  data,  5,y  are  the  components  of  the  deviatoric  stress  tensor,  and  ar 
is  the  material’s  yield  stress.  Strain  hardening  and  softening  is  simulated  by 
making  the  yield  stress  a  function  of  plastic  strain.  In  the  uniaxial  case,  5,-y  = 


2 


and  Equation  (14)  becomes 


Since  the  overstress  law  only  calculates  viscoplasticity  when  the  applied 
stress  exceeds  the  yield  stress,  another  law  is  needed  to  account  for  creep  at 
lower  stress  levels.  Norton’s  Law  for  Secondary  Creep  is  a  special  case  of  the 
Power  Law  in  which  a  constant  internal  state  is  assumed,  making  creep  strain 
rate  dependent  only  upon  applied  stress  with  no  strain  hardening  or  softening. 
Creep  strain  rate  is  predicted  to  be 


where  7C  and  $  are  material  constants  determined  from  creep  tests  conducted 
at  two  different  stress  levels,  and  ff0  is  a  normalizing  stress  generally  given  an 
arbitrary  value  such  as  100  KSI.  In  the  uniaxial  case,  Equation  (16)  becomes 


:=7 

*0 


The  numerical  Euler  time  integration  is  less  involved  than  that  required  for 
the  Bodner-Partom  model,  consisting  of  the  following  steps  performed  during 
eac  h  time  step: 


The  material  parameters  used  in  this  investigation  for  the 
Ove  rs  tre  ss/ Nor  ton ’s  Law  model  were  calculated  from  cyclic  test  data  and  are 


shown  in  Table  II. 


Table  II.  Overstress  /  Norton’s  Law  Model  Material  Constants 

Material:  Inconel  718  @1200  F 

K 

0.0107 /sec 

n 

1.0 

I* 

1.520xl0"2Vsec 

*0 

1.0  PSI 

0 

4.022 

Time  Step  Selection 


For  stable,  accurate  results,  the  time  increment  dt  must  not  exceed  the 
maximum  allowable  value  for  the  model  used.  Cormeau  (Ref  13)  shows  how  to 
determine  maximum  time  step  sites  for  Euler  time  integration  schemes  by 
using  matrix  algebra  to  find  eigenvalues  for  the  simultaneous  differential  equa¬ 
tions  involved.  In  two  or  more  dimensions  the  calculations  are  somewhat 
cumbersome,  however  in  the  uniaxial  case  the  procedure  simplifies  to 

(21) 

where 


X  =  £7/ 


and 


H-. 


di‘ 


da 


For  the  uniaxial  Overstress  and  Norton  models. 


(22) 

(23) 


dt. 


dt , 


2(7  Y 
lE 


( Ove  rstress,  n  =1) 


(24) 


I'EHafr') 


(  No  r  t  o  n  ) 


(25) 


Maximum  time  steps  for  the  Bodne  r-Partom  model  can  also  be  calculated 


for  a  given  stress  level: 


(26) 


In  general,  the  maximum  allowable  time  step  during  viscoplastic  deforma¬ 
tion  will  be  several  orders  of  magnitude  smaller  than  the  maximum  allowable 
time  step  during  creep.  For  stress-strain  predictions  in  which  stresses  or  strain 
rates  vary  between  high  and  low  values  it  is  therefore  often  desireable  to  vary 
the  time  step  size  during  program  execution  in  order  to  reduce  the  total 
number  of  iterations  required.  Hinnerichs  (Ref  10)  describes  a  method  for  vary¬ 
ing  the  time  step  size  based  on  the  changes  in  stress  and  strain  which  occur 
during  a  given  time  step.  For  simplicity,  FINELS  does  not  vary  the  time  step, 
although  such  a  method  could  be  incorporated  if  necessary. 


Additional  Stability  Oiteria 

For  integration  of  a  model  which  contains  an  evolving  internal  state  vari¬ 
able,  proper  time  step  selection  is  necessary  but  not  sufficient  for  ensuring 
numerical  stability.  Ponter  (Ref  14)  shows  that  for  stability, 

d<r  d  i  p  —  dZ  dZ  >0  (27) 

during  each  time  step,  where  do,  dip ,  dZ ,  and  dZ  represent  increments  in 
stress,  inelastic  strain  rate,  internal  state  variable,  and  evolution  rate  of  the 
internal  state  variable,  res  pe  c  t  ive  I  y . 


In  general,  viscoplas  tic  lty  models  will  be  unconditionally  stable  when  strain 
hardening  is  modelled  (Refs  7,14).  Such  models  will  only  become  conditionally 
stable  when  recovery  rates  are  large.  Strain  softening  models,  on  the  other 
hand,  are  only  conditionally  stable  when  the  increments  in  stress  and  strain  rate 
are  large  compared  with  those  of  the  internal  variable  and  its  evolution  rate 
(Refs  7,14).  This  limits  the  ability  of  viscoplasticity  models  to  predict  the 
behavior  of  strain  softening  materials,  especially  at  high  applied  stresses  and 
low  frequencies  where  the  model  is  most  likely  to  become  unstable.  This  did  in 
fact  occur  and  is  documented  in  Section  IV. 
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IV.  Results 


Program  Validation 

The  program  used  for  this  investigation  was  FINELS  (See  Appendix  C),  a 
FORTRAN  code  developed  to  model  viscoplasticity  and  creep  in  systems  of  uniax¬ 
ial  bars  using  finite  element  methods. 

Atime  step  convergence  study  was  undertaken  for  FINELS  in  order  to  verify 
the  program’s  accuracy  and  to  validate  the  maximum  time  step  size  require¬ 
ment.  The  test  problem  described  by  Mercer  (Ref  9)  was  chosen,  in  which  the 
program  predicts  the  behavior  of  a  uniaxial  bar  made  of  IN  100  subjected  to 
cyclic  stress  of  ±200  KSI  peak  values  applied  at  a  frequency  of  10  cpm.  The 
problem  was  run  using  time  steps  of  0.1,  0.05,  0.01,  0.005,  and  0.001  seconds. 
For  a  time  step  dt  of  0.1  seconds,  pronounced  ratcheting  was  evident  as  the 
stress-strain  loop  shifted  in  the  direction  of  decreasing  strain  with  each  applied 
stress  cycle.  For  dt  =  0.05  seconds,  the  ratcheting  was  no  longer  evident 
although  the  stress-strain  curve  does  not  match  Mercer’s.  For  dt  =  0.01 
seconds,  the  curves  match  more  closely  and  for  dt  =  0.005  seconds  and  dt  = 
0.001  seconds,  the  FINELS  curve  is  almost  identical  to  Mercer’s  (See  Figures  1 
through  5). 

The  ratcheting  phenomenon  appears  to  be  related  to  the  numerical  stabil¬ 
ity  of  the  model.  Using  the  Bodner-Partom  material  constants  shown  in  Table  I 
for  IN  100,  for  ±200  KSI  peak  stress  the  maximum  allowable  time  step  given  by 
Equation  26  is 

|Umax=0.009  seconds 

The  inaccuracies  and  ratcheting  shown  for  dt  =0  01  seconds,  0.05  seconds,  and 
0.1  seconds  respectively  indicate  that  the  model  is  not  stable  for  these  time 
step  sizes.  The  maximum  tune  step  sue  therefore  appears  to  be  valid  since 


$ 

.*>«* 


:u 

1 


\\  -o  £ 


m 

O 

m 

o 

o 

o 

o 

o 

o 

+ 

"+■ 

■+• 

LJ 

LU 

Ld 

O 

o 

O 

• 

CM 

6 

CN 

(iSd)  SS3diS 


v*y!v!v* 


w  Q. 
o 

*?  «53 

hu 

nnZ 

•o 

c  —  i 
o  ®  o 
*o  £ 
0)  o  co 
®  E 

c  c  w  i 


■D  « 


CO  2 
0)  C 
C  O 


<u  O  o 

*E  *o  ."2  ’'S 

W  us  —  <n 
OoOu 
Q  > 


O  C 

co  a. 


o 

o 

CM 

ii 

o\ 
o  + 


II  CO 
CO 
>>  (U 
O  w. 
c  +~> 

s « 

H 


« s 

ia 


wWw.iw&iWwwft 


ww.*' 


2.0E+005 


* 

o 

JC 

CO 

w 


E 

o 

k_ 

o» 

o 

L. 

Q. 


Figure  4.  IN  100  Time  Convergence  Study,  DT  =  0.005  sec. 


inaccuracies  or  ratcheting  did  not  occur  when  the  time  step  used  was  smaller 
than  0.009  seconds. 

Asimilar  time  convergence  study  was  performed  for  Inconel  718  at  1200  F 
subjected  to  uniaxial  cyclic  strains  of  ±0.008  peak  values  applied  at  a  frequency 
of  10  cpm.  In  this  case,  cyclic  strain  was  analyzed  due  to  the  model’s  inade¬ 
quate  cyclic  stress  prediction.  The  time  steps  used  were  0.1  seconds,  0.05 
seconds,  0.01  seconds,  and  0.005  seconds.  The  results  for  each  time  step  size 
are  shown  in  Figures  6  through  8. 

Under  cyclic  strain,  no  ratcheting  is  produced.  However,  convergence  was 
noted  since  the  stress-strain  curve  for  dt  =0.005  seconds  was  closely  matched 
by  that  for  dt  =0.01  seconds  and  less  closely  matched  as  the  time  step  sizes 
inc  re  ased . 

The  maximum  allowable  time  step  size  can  be  calculated  using  Equation  26. 
For  a  peak  stress  of  approximately  150  KSI  it  is 

dtmtx  =  0.006  seconds 

Therefore  steadily  increasing  divergence  is  expected  for  dt  =0.01  seconds,  0.05 
seconds,  and  0.1  seconds,  respectively. 

For  this  investigation,  Cormeau’s  method  for  determining  the  maximum 
allowable  time  step  for  numerical  stability  during  Euler  time  integration  was 
found  to  be  valid  for  the  Bod  n  e  r-Parto  m  model  as  well  as  for  the  other  viscoplas¬ 
ticity  models  used.  However,  if  conditions  other  than  the  ones  used  in  this 
investigation  are  to  be  modelled,  a  limited  time  convergence  study  may  be 
desireable  to  confirm  the  accuracy  of  the  model 

Accuracy  is  also  affected  by  the  number  of  time  steps  per  cycle  of  applied 
stress  or  strain.  No  systematic  attempt  was  made  to  calculate  the  minimum 
number  of  time  steps  per  cycle  necessary  for  accurate  results.  However,  at 
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Figure  7.  IN  718  Time  Convergence  Study,  DT  =  0.05  sec. 


high  stress,  the  maximum  time  steps  allowed  by  Equations  24  through  26  are 
very  small  compared  to  the  cyclic  period.  This  placed  a  lower  limit  upon  the 
number  of  time  steps  which  occurred  during  each  cycle.  It  was  found  that  this 
resulted  in  a  sufficient  number  of  time  steps  per  cycle  to  provide  accurate 
results.  The  time  step  sizes  used  for  predictions  at  the  highest  stress  levels 
were  therefore  used  for  lower  stress  levels  as  well. 

FINELS  was  also  used  to  duplicate  the  results  obtained  by  Hinnerichs  and 
Palazotto  (Ref  15)  for  a  three  bar  linkage  undergoing  viscoplasticity  and  creep. 
This  is  discussed  in  Appendix  B. 

Comparison  of  FINELS  Output  with  Test  Data 

FINELS’  accuracy  was  further  confirmed  by  duplicating  the  results  Stouffer 
(Ref  11)  obtained  when  determining  the  Bodner-Partom  material  constants  for 
IN  100  at  1350  F  from  uniaxial  tensile  and  creep  test  data  (See  Appendix  A). 
FINELS’  results  matched  Stouffer’s;  the  stresses  and  strains  calculated  using 
the  Bodner-Partom  model  were  in  close  agreement  with  tensile  test  data,  but 
showed  the  same  discrepancies  as  Stouffer’s  results  when  compared  to  creep 
test  data. 

Using  Stouffer’s  constants  for  IN  100,  the  Bodner-Partom  model  calculates 
creep  strain  rates  which  at  low  stress  levels  are  initially  too  low,  and  which  at 
higher  stress  levels  are  initially  too  high,  although  in  all  cases  the  creep  strain 
rate  eventually  stabilizes  at  the  correct  value.  This  behavior  is  directly  related 
to  the  thermal  recovery  of  hardening  defined  in  Equation  8,  which  controls  the 
rate  at  which  the  parameter  Z  evolves.  In  this  case,  Z  does  not  evolve  rapidly 
enough  to  allow  creep  strain  rate  to  stabilize  quickly.  If  Z  is  initially  too  high, 
the  creep  rate  will  initially  be  too  low  until  Z  evolves  to  its  proper  value.  Simi¬ 
larly,  if  Z  is  initially  too  low.  the  creep  rate  will  initially  he  too  high.  This  will 
lead  to  cumulative  error  over  time  if  the  rate  at  which  Z  evolves  is  too  slow. 


Stouffer  points  out  that  the  material  constants  related  to  thermal  recovery, 
A,  r,  and  Zo ,  are  not  well-defined  for  IN  100  due  to  the  large  scatter  in  the  creep 
test  data.  This  behavior  is  therefore  a  potential  source  for  cumulative  error 
when  modelling  IN  100  in  situations  where  significant  creep  will  occur. 

Computer  Accuracy.  It  was  initially  thought  that  the  accuracy  of  the  computer 
may  have  caused  the  inaccuracy  for  low  stress  levels.  The  VAX  1  1/780  upon 
which  FINELS  was  run  reads  any  number  less  than  1.0X10-38  as  zero.  Since  ine¬ 
lastic  strain  rate  contributes  to  the  evolution  rate  of  Z,  and  since  the  initial  ine¬ 
lastic  strain  rate  is  small  for  low  stress,  it  was  thought  that  the  computer  may 
have  been  causing  errors  by  reading  the  strain  rate  contribution  as  zero. 

Confirmation  was  attempted  by  modifying  FINELS  to  change  the  value  of  the 
number  read  by  the  computer  as  zero.  It  was  thought  that  artificially  altering 
the  computer’s  accuracy  in  this  manner  would  further  slow  the  rate  at  which  Z 
evolves,  thereby  increasing  the  time  required  for  the  creep  rate  to  stabilize  at 
its  prope  r  value  . 

It  was  found  that  this  was  not  the  primary  cause  of  the  slow  creep  rate  evo¬ 
lution,  since  an  increase  of  the  "zero"  value  from  1  OxlO~38  to  1.0X10-26  did  not 
cause  significant  differences  in  the  predictions  (See  Figure  9). 

It  was  also  found  that  further  increases  in  "zero"  value  to  1.0X10-24, 
1.0X10-23,  and  l.OxlO-20  did  significantly  affect  the  accuracy  of  the  predictions, 
however  (See  Figure  9).  This  is  not  due  to  inaccuracies  in  the  evolution  rates  of 
the  internal  variables.  Rather,  it  is  because  the  exponential  term  used  in  Equa¬ 
tion  2  to  calculate  the  viscoplastic  strain  rate  has  a  value  which  is  very  small  (on 
the  order  of  1.0X10-'3)  n  n  <  1  »■  r  these  conditions  even  when  the  internal  state  vari¬ 
able  reaches  its  steady-state  value.  This  represents  a  source  of  progressive 
error  when  performing  creep  predictions  using  computer  systems  whose  accu¬ 
racy  is  less  than  that  of  the  V  \X  |  |  7 M) 


Conditions:  Uniaxial  tension 
(Constant  stress) 


Peak  Stress  (KS1) 


Table  III 

Matrix  of  Predicted  Results 
Bodne  r-Partom  Model 
Un iaxial  Cy c  lie  Stress 
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gure  10.  IN  718  Prediction  (Table  I  Material  Properties,  Initial  Set) 


Bodner  model:  Overstress  model: 


m 

O 

m 

o 

O 

o 

o 

o 

o 

■+■ 

4* 

+ 

LU 

LxJ 

LxJ 

O 

O 

O 

• 

CM 

• 

o 

CM 

1 

(iSd)  SS3HJ.S 


o 

O  O 

TJ  TJ 
O  —  O  <U 
"o  o  a>  u 
‘_P  -w  Q 
co 

CO  V_  — 

S'"  II  ® 
II  o 

>  CO  c 
O  w  g>  b 

£i!o> 

"eg.8 

“OW  m 
o 


Figure  15.  Uniaxial  Cyclic  Strain  Prediction  (0.01  CPM) 


The  revised  Bodne  r-Partom  material  constants  were  calculated  from 
Inconel  718  cyclic  test  data  using  the  methods  described  by  StoufFer  (See  Table 
I).  When  predicting  strain-control  behavior,  they  were  found  to  portray  the 
material’s  viscoplastic  behavior  reasonably  well  over  the  range  of  frequencies 
modelled  (See  Figs.  11  through  16).  However,  as  with  the  original  material  con¬ 
stants,  they  produced  stress-strain  curves  with  sharply  defined  yielding  points, 
rather  than  the  smooth  stress-strain  curves  produced  by  the  test  data.  This  is 
due  to  the  exponential  relationship  between  stress  and  viscoplastic  strain  rate 
assumed  in  the  Bodne  r-Partom  model  (See  Equations  1  through  4). 

Using  the  new  constants,  the  Bodner-Partom  model  did  not  predict  stress- 
control  behavior  adequately.  No  appreciable  viscoplasticity  was  predicted 
although  the  test  data  indicate  that  considerable  viscoplasticity  occurs  under 
the  conditions  modelled  (See  Figures  17  through  19). 

It  was  also  found  that  when  strain  softening  behavior  was  closely  modelled, 
pure  creep  was  not  adequately  portrayed.  This  is  because  the  model  does  not 
allow  Z  to  drop  to  a  low  enough  value  for  significant  creep  to  occur;  in  Equation 
5,  when  Z0  is  greater  than  Z Zx  becomes  the  lower  limit  for  Z.  Thermal 
recovery  of  plastic  work  cannot  cause  Z  to  become  less  than  Z\  in  this  case.  It 
is  possible  that  this  could  be  corrected  by  proper  selection  of  recovery  con¬ 
stants;  however,  this  was  not  attempted. 

Under  stress  control,  the  model  became  unstable  at  high  applied  stresses 
and  low  frequencies.  Instabilities  occured  at  frequencies  of  10  cpm  and  lower 
for  peak  stresses  of  ±150  KSI  and  at  frequencies  of  0.1  cpm  and  lower  for  peak 
stresses  of  ±120  KSI.  As  applied  peak  stresses  decreased,  the  range  of  frequen¬ 
cies  in  which  instability  occurred  grew  smaller  This  is  <  onsistant  with  the  sta¬ 
bility  requirements  described  by  Equation  27  At  high  stress  and  low  frequen¬ 
cies,  the  increments  in  state  variable  and  its  evolution  rate  during  a  given  time 


Figure  18.  Uniaxial  Cyclic  Stress  Prediction  (150  KSI,  100  CPM) 
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step  become  large,  and  for  a  strain  softening  material  are  sufficent  to  violate 
the  requirements  for  stability 

Table  III  also  shows  the  conditions  for  which  the  various  types  of  Bodner- 
Partom  model  rest '>1.  es  occurred 

Frequency  Response  of  the  Overstress/ Norton  Model 

The  same  conditions  were  modelled  {See  Table  IV)  using  a  version  of  FINELS 
modified  to  use  a  combination  of  Overstress  law  to  determine  viscoplasticity  and 
Norton’s  Law  for  Secondary  Creep.  Strain  softening  was  modelled  by  making  the 
yield  stress  in  the  Overstress  equation  a  function  of  accumulated  plastic  strain: 

o  > ■  =  1 00  KS  /  {0.9+0.2expj  — mVV,,])  (28) 

where  YV^?  is  the  net  plastic  work  and  is  given  by 


or,  during  time  step  i. 
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Table  IV 

Matrix  of  Predicted  Results 
Ove  rs  tress  /  Nor  ton  Model 
Uniaxial  Cvc  lie  Stress 
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The  Overstress  Norton  model  also  yielded  good  predictions  for  strain- 
control  behavior  (See  Figures  11  through  16).  In  fact,  for  high  frequencies  (1 
cpm  and  higher),  the  Overstress  /  Norton  model  portrayed  the  material’s 
behavior  slightly  more  accurately  than  the  Bodner-Partom  model,  since  the 
Overstress/ Norton  model  allows  a  smoother  transition  between  creep  and  visco¬ 
plasticity  at  high  stress  <  strain  application  rates.  However,  at  lower  frequencies 
the  "on-off"  nature  of  the  Overstress  law  causes  the  stress-strain  curve  to 
assume  the  more  squared-off  shape  characteristic  of  the  Bodner-Partom  model. 
This  is  because  when  the  stress  reaches  the  yield  stress,  the  Overstress  law 
causes  the  inelastic  strain  rates  to  increase  by  several  orders  of  magnitude  over 
a  relatively  small  increment  in  stress. 

At  the  highest  frequency  modelled,  100  cpm,  a  discrepancy  was  noted 
where  the  Ove  rstre  ss  /  Nor  ton  model  predicted  considerably  less  plastic  strain 
than  the  Bodner-Partom  model.  In  this  case,  although  the  Overstress/Norton 
model  begins  yielding  at  a  lower  stress,  the  Bodner-Partom  model  produces  a 
higher  viscoplastic  strain  rate  once  viscoplasticity  begins. 

The  Overstress/ Norton  model  did  not  adequately  portray  stress-control 
behavior  (See  Figures  17  through  19).  Its  results  were  similar  to  those  of  the 
Bodner-Partom  model;  although  some  viscoplasticity  was  shown  at  120  KSI,  the 
highest  applied  stress  for  which  test  data  were  available  (See  Figure  17).  the 
predictions  did  not  agree  with  the  test  data,  while  no  appreciable  viscoplasticity 
was  shown  in  any  other  case  (See  Figure  18).  At  low  stress,  however,  Norton's 
Law  provided  good  predictions  of  creep  response  (See  Figure  19). 

It  is  possible  that  the  Overstress  iaw  could  be  made  to  more  closely  match 
the  material  behavior  bv  making  the  yield  stress  a  function  of  total  strain,  in 
addition  to  plastic  work,  and  by  including  recovery  terms  in  a  manner  similar  to 
that  of  the  Bodner-Partom  method.  Th  is  is  an  area  which  requires  further  studv 


In  addition,  the  Overstress  model  was  found  to  be  unstable  at  high  applied 


peak  stress  and  low  frequency  during  stress  control.  Instability  occurred  under 


the  same  conditions  as  for  the  Bodner-Partom  model.  This  is  also  consistant 


with  the  stability  criteria  described  by  Equation  27.  In  this  case,  the  yield  stress 


evolves  in  a  manner  similar  to  that  of  the  state  variable  Z  in  the  Bodner-Partom 


model.  It  is  therefore  to  be  expected  that  instability  should  occur  under 


approximately  the  same  conditions  (high  stress,  low  frequency)  as  for  the 


Bodner-Partom  predictions. 


Table  IV  also  shows  the  conditions  at  which  the  various  types  of  behavior 


occurred,  i.e.,  instability,  viscoplasticity,  creep,  and  primarily  elastic  response 


Here  it  should  be  noted  that  the  occurence  of  instability  indicates  that 


viscoplasticity  is  occurring,  since  instability  is  triggered  by  evolution  of  the 


internal  state  variable  (as  shown  in  Equation  27),  which  in  turn  is  brought  about 


by  the  accumulation  of  viscoplastic  strain. 


C(yclic  Oeep  (IN  100) 


It  has  been  noted  that  the  isotropic  strain-softening  Bodner-Partom  model 


is  unable  to  model  creep.  This  does  not  mean  that  the  Bodner-Partom  model 


cannot  model  creep  under  any  conditions.  In  order  to  confirm  that  the  Bodner- 


Partom  model  can  model  cyclic  creep,  a  typical  uniaxial  cyclic  creep  response 


(100  KSIat  0.001  cpm)  was  modelled  using  the  constants  for  IN  100  (See  Table  I) 


Since  IN  100  is  strain-hardening,  the  condition  which  prohibits  calculation  of 


creep  in  a  s  tr  ain-soft  e  n  in  g  material  does  not  occur. 


The  Bodner-Partom  model  results  were  compared  to  those  of  the 


Overstress  Norton  model  for  the  same  conditions.  In  this  case  the  material 


stants  used  were  (Ref  10): 
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Stress  Concentration 


borne  by  the  remaining  two  unyielded  members.  If  the  applied  load  is  increased 
sufficiently,  a  second  member  will  begin  yielding,  and  finally  the  third.  The 
model  therefore  shows  how  the  inelastic  behavior  of  one  element  affects  the 


behavior  of  adjacent  elements.  In  this  manner,  the  model  simulates  the  forma¬ 
tion  and  growth  of  a  plastic  zone  around  a  crack  tip  or  other  area  of  stress  con¬ 
centration  in  a  two-dimensional  model. 

Aplot  of  typical  three-bar  linkage  behavior  over  time  during  cyclic  applied 
stress  is  shown  in  Figure  22.  In  it,  the  points  where  load  is  transferred  from  one 
element  to  others  due  to  viscoplastic  response  can  clearly  be  seen 

In  general,  the  behavior  of  the  yielding  elements  (elements  #2  and  #3)  was 

neither  that  of  pure  stress  control  nor  that  of  pure  strain  control  (see  Figures 

23  through  25).  However,  for  applied  forces  of  200  KPS  and  250  KPS  peak  value, 
the  elements'  behavior  was  that  of  pure  stress  control  at  the  highest  and  lowest 
frequencies,  where  the  peak  stresses  remained  the  same  while  the  peak  strains 
grew  greater  during  each  successive  cycle.  At  intermediate  frequencies,  typi¬ 
cally  10  cpm  and  1  cpm,  some  behavior  reminiscent  of  strain  control  was  noted, 
i.e.  high  initial  peak  stress  and  lower  peak  stresses  during  subsequent  loadings 
(See  Figure  24). 

The  Bodner-Partom  predictions  were  then  compared  to  predictions 

obtained  for  the  same  conditions  using  the  Ove  rs  tre  ss  /  Nor  ton  model  (See  Table 
VI).  The  Ove  rs  tress  /  Nor  ton  model  predicts  behavior  similar  to  that  of  the 

Bodner-Partom  model  although  the  resulte  do  not  always  match  well.  A  com¬ 
parison  of  results  is  shown  in  Figures  23  through  25. 
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Figure  23.  Stress  Concentrotion  Predictions  (250  KPS,  100  CPM) 
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figure  25.  Stress  Concentration  Predictions  (100  KPS,  0.001  CPM) 


Table  VI 

Matrix  of  Predicted  Results 
Ove  rs  tress/ Nor  ton  Model 


Sim  ulate  d  Stress  C 

lone  e  ntrati 

r\  1  t— t  t  r/-n  r<  \ 

Freauencv  l 

Peak  Force  (KPS) 

100 

10 

1 

0.1 

250 

E 

\P 

VP 

VP 

200 

E 

VP 

VP 

VP 

150 

E 

E 

E 

E 

100 

E 

E 

E 

E 

KEY. 

I  =Instability 

VP  =Viscoplastic  Response 
C=  Creep  Response 

E  =Elastic  Response  (Little  or  no  vise  op  las  tic  ity  or  creel 


The  Ove  rs  tre  ss/ Nor  ton  model  predicted  lower  peak  strains  than  those  of 
the  Bodner-Partom  model  for  an  applied  force  of  250  KPS  (See  Figure  23)  while 
predicting  higher  peak  strains  for  the  corresponding  frequencies  when  the  force 
was  200  KPS  (See  Figure  24).  In  addition,  at  low  applied  forces  and  frequencies, 
the  models’  predictions  diverged  due  to  the  Ove rstress/ Norton ’s  Law  model’s 
portrayal  of  creep  and  the  Bodner-Partom  model’s  previously  noted  inability  to 
portray  creep  under  these  conditions  (See  Figure  25). 

Instability  did  not  occur.  However,  it  must  be  noted  that  conditions  of  high 
applied  stress  and  low  frequency  produce  results  which  cannot  be  considered 
approximations  of  real-life  results,  since  the  strains  produced  by  maintaining 
high  loads  for  long  periods  are  well  past  the  rupture  limit  for  the  actual 
material.  Also,  the  assumption  of  small  deformations  no  longer  applies.  These 
predictions  are  therefore  of  questionable  value  for  comparing  the  models' 
behavior. 

Tables  V  and  VI  also  show  the  conditions  for  which  the  various  types  of 
behavior  occurred,  i.e  .  .vise  op  las  tic  tty ,  creep,  and  primarily  elastic  response 

It  is  perhaps  more  informative  to  view  the  results  in  terms  of  ac  c  u  m  it  la  t  e  <1 
inelastic  strain  vs.  frequency  as  shown  in  Figure  26.  This  better  shows  the 


conditions  for  which  a  viscoplasticity  model  is  necessary  and  those  for  which  a 
simpler  creep  model  or  even  a  purely  elastic  model  will  suffice.  At  high  stresses, 
inelastic  strains  are  largely  due  to  viscoplasticity.  At  low  stresses,  inelastic 
strains  are  largely  due  to  creep.  Where  little  or  no  inelastic  strains  occur,  the 
response  is  primarily  elastic. 


V'.  Conclusions 


Summary  of  W>rk  Done 

The  cyclic  behavior  of  the  Bodne  r-Partom  and  combined 
Overstress/ Norton ’s  Law  models  were  investigated  over  a  range  of  frequencies 
and  applied  cyclic  stresses  and  strains  using  the  viscoplasticity  prediction  code 
F1NELS  developed  for  the  purpose.  For  simplicity,  only  isotropic  hardening  was 
investigated,  although  kinematic  hardening  can  also  be  modelled.  F1NELS  was 
validated  by  comparison  of  its  results  with  those  of  other  codes  and  with  cyclic 
and  non-cyclic  tensile  and  creep  test  data.  Its  ability  to  simulate  stress  concen¬ 
tration  was  tested  by  running  the  three-bar  linkage  problems  formulated  by  Hin- 
nerichs  and  Palazotto  (Ref  15).  FINELS  was  used  to  check  the  validity  of 
Cormeau’s  (Ref  13)  requirement  for  the  maximum  allowable  time  step  size  for  a 
stable  model.  Uniaxial  cyclic  behavior  was  predicted  for  cyclic  applied  stresses 
and  strains  comprising  a  range  of  frequencies  and  amplitudes.  Cyclic  stress 
concentration  behavior  was  predicted  by  mode  Ring  a  three-bar  linkage  sub¬ 
jected  to  cyclic  applied  loads  comprising  a  range  of  frequencies  and  amplitudes. 

Characteristics  of  the  Bodner-Partom  Model 

Frequency  Response.  A  viscoplasticity  model  is  required  to  predict  material 
behavior  at  high  stresses  (alternately,  high  strain  rates)  over  a  range  of  frequen¬ 
cies,  and  at  intermediate  stresses  (or  strain  rates)  at  low  frequencies  ( i . e . ,  con¬ 
ditions  corresponding  to  the  upper  right  region  of  Tables  III  through  VI),  since 
these  are  the  conditions  near  a  crack  tip  for  which,  in  real  life,  large  viscoplastic 
strains  will  occur. 

Both  the  Bodner-Partom  model  and  the  Overstress  model  can  be  used  to 
predict  viscoplastic  response  although  their  predicted  results  diverge  somewhat 
at  the  highest  applied  stresses  and  strain  rates. 


Viscoplasticity  models  are  not  always  required,  however.  For  low  stresses 
applied  at  low  frequencies  (the  lower  right  region  of  Tables  III  through  VI),  the 
material’s  behavior  is  largely  due  to  creep.  In  these  cases,  a  simple  model  such 
as  Norton’s  Law  will  adequately  predict  material  behavior.  For  a  strain-softening 
material,  Norton’s  Law  is  actually  preferable  to  the  Bodne r-Partom  model  due  to 
the  inability  to  predict  creep  using  the  isotropic  strain-softening  Bodner-Partom 
model  under  these  conditions. 

Also,  for  intermediate  and  low  stresses  at  high  frequencies  (the  lower  left 
region  of  Tables  III  through  VI),  little  or  no  viscoplasticity  or  creep  will  occur. 
Under  these  conditions,  a  purely  elastic  model  will  suffice  for  stress-strain  pred- 
ic  tion . 

Greep  Behavior.  Inadequate  creep  test  data  prevents  calculation  of  accurate 
recovery  parameters.  This  in  turn  prevents  accurate  modelling  of  the  evolution 
rate  for  Z.  Although  Z  will  eventually  evolve  to  its  appropriate  value,  a  poten¬ 
tially  large  initial  error  will  be  introduced  when  creep  predictions  are  made. 

Kinematic  Hardening.  It  appears  that  Inconel  718  shows  significant  kinematic 
hardening  as  well  as  isotropic  hardening.  Beaman  (Ref  8)  showed  how  the 
assumption  of  kinematic  hardening  produces  results  which  match  cyclic  test 
data  much  more  closely  than  results  produced  by  purely  isotropically  hardening 
models.  Lindholm,  et  al  (Ref  7)  also  showed  how  the  addition  of  kinematic  har¬ 
dening  terms  to  Z  allows  the  Bodner-Partom  model  to  predict  smoother,  more 
accurate  e  las  tic -vise  op  las  tic  transition. 

Time  Step  Siie.  Ratcheting  and  other  cumulative  errors  can  be  eliminated  by 
proper  selection  of  the  time  step  size.  The  time  step  size  is  critical  to  ensuring 
numerical  stability  of  the  Bodner-Partom  model.  Cormeau’s  approach  to  deter¬ 
mining  the  maximum  allowable  time  step  was  found  to  apply  to  the  Bodner- 
Partom  model  as  well  as  to  others  When  the  time  step  exceeds  this  maximum 
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allowable  value,  cumulative  error  is  introduced.  Further  increases  in  time  step 


D 


& 


size  cause  ratcheting  and  eventually  cause  the  model  to  "blow  up”  numerically. 

Strain  Softening.  It  has  been  previously  noted  that  the  Bodne r-Partom  model 
does  not  adequately  portray  creep  when  a  strain-softening  material  such  as 
Inconel  718  is  modelled  assuming  isotropic  hardening  only.  In  Equation  5,  when 
Z o  is  greater  than  Z\,  Zx  becomes  the  lower  limit  for  Z.  Thermal  recovery  of 
plastic  work  cannot  cause  Z  to  become  less  than  Zx  in  this  case.  The  Bodner- 
Partom  model  can  portray  creep  adequatetely  when  isotropic  strain-hardening 
materials  such  as  Inconel  100  are  modelled. 

Numerical  Instability.  Both  models  share  a  common  limitation  with  all  strain 
softening  models:  that  of  conditional  stability.  In  this  case,  however,  it  is  not  of 
practical  consequence  since  the  conditions  under  which  it  occurs  (high  stresses 
and  low  frequencies;  in  this  case,  10  cpm  and  lower  at  150  KSI,  1  cpm  and  lower 
at  120  KSI)  are  those  under  which,  in  real  life,  the  material  would  rupture. 


Simulated  Stress  Concentration. 


The  behavior  of  Inconel  718  at  a  stress  concentration  when  subject  to  cyclic 
remote  stress  or  displacement  is  generally  closer  to  that  of  pure  stress  control 
than  to  that  of  pure  strain  control,  although  under  conditions  of  intermediate 
stress  and  frequency  behavior  similar  to  strain  control  can  occur,  namely  high 
stresses  during  the  first  cycle  and  lower  stresses  during  subsequent  cycles. 


*  M  *.»“**  ' 


VI.  Recommendations 


Oeep  in  a  Strain  Softening  Material. 

The  isotropic  Bodner-Partom  formulation  for  thermal  recovery  does  not 
allow  creep  rate  in  a  strain  softening  material.  Alternate  methods  should  be 
developed  for  determining  the  evolution  of  the  state  variable  Z  to  allow  it  to 
evolve  below  the  lower  limit  imposed  upon  it  by  Equation  (5)  for  strain  softening 
m  aterials. 

Kinematic  Hardening. 

Kinematic  hardening  should  be  accounted  for  in  all  cyclic  viscoplastic  ana¬ 
lyses  involving  high-temperature  turbine  materials. 

Continued  Investigation. 

The  next  step  in  an  ongoing  investigation  should  be  to  extend  analysis  to 
problems  in  two  and  three  dimensions.  ITie  results  of  this  investigation  could  be 
used  to  help  identify  the  conditions  at  a  crack  tip  for  which  each  type  of 
behavior  (viscoplasticity,  creep,  elasticity)  can  be  expected  to  occur. 
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Appendix  A 

Comparison  of  FINE!  3  Results  with  Stouffer’s  Predictions 


Stouffer  (Ref  11)  showed  how  to  obtain  Bod  n  e  r-Par  to  m  constants  for  a 
material  using  data  obtained  from  uniaxial  tensile  and  creep  tests.  The  con¬ 
stants  thus  obtained  were  then  validated  by  comparing  Bodner-Partom  predic¬ 
tions  with  the  test  data  used  to  obtain  the  constants.  FINELS  was  run  using 
Stouffer’s  constants  and  the  results  compared  to  Stouffer’s  results  and  test 
data . 

Stouffer  calculated  the  Bodner-Partom  constants  for  IN  100  to  be  those 
shown  in  Table  I.  The  tensile  test  data  was  obtained  from  uniaxial  tensile  tests 
conducted  using  applied  strain  rates  ranging  from  1 .4  2  Xl  0-3/' s  e  c  to 
1 .67XI0-8/ sec  .  The  creep  test  data  was  obtained  from  uniaxial  tests  conducted 
using  applied  stresses  ranging  from  72  KSI  to  130  KSI.  Stouffer  predicted 
responses  for  applied  stresses  of  70  KSI,  90  KSI.  110  KSI,  120  KSI,  and  130  KSI, 
and  for  applied  strain  rates  of  1 .4X10_3/sec  ,  5 .5  Xl  0“6/ se  c  ,  and  6  .6X10~®  /  se  c  . 
The  material  behavior  was  predicted  again  using  FINELS  with  the  same  applied 
stresses  and  strain  rates,  and  compared  with  Stouffer’s  predictions  and  test 
data  as  shown  in  Figures  A-l  and  A-2. 


Comparison  of  FIN  ELS 
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Appendix  B 

The  Three -Bar  Linkage 

The  three-bar  linkage,  as  originated  by  Kawahara  (Ref  16)  and  expanded 
upon  by  Hinnerichs  and  Palazotto  (Ref  15),  is  a  method  for  simulating  material 
behavior  at  a  stress  concentration  by  using  uniaxial  bar  elements. 

The  three-bar  linkage  consists  of  three  uniaxial  bar  elements  as  shown  in 
Figure  21.  Stress  concentration  is  simulated  by  applying  a  load  F  at  the  central 
node,  oriented  such  that  the  force  acting  upon  one  element  is  greater  than 
those  acting  upon  the  other  two.  Viscoplasticity  and  creep  can  be  included  as 
described  elsewhere  in  this  paper. 

FINELS  was  used  to  duplicate  the  results  obtained  by  Hinnerichs  and  Pala¬ 
zotto  (Ref  15)  when  they  modelled  the  effects  of  viscoplasticity  and  creep  in  a 
three-bar  linkage.  The  constants,  time  increments,  and  applied  forces  were  cal¬ 
culated  in  the  manner  shown  in  Ref  15  and  are  shown  in  Table  VII.  The 
material’s  static  stress-strain  slope  in  the  plastic  region  was  modelled  by  assum¬ 
ing  the  yield  stress  to  be  linearly  related  to  total  strain,  in  this  case 

ar  =  [(Tr]o  +  0.1  E(  |c  |  —  |c  y  |o)  (  le  l>l£  y  lo  ) 

=  l*7  y]  o  (  lf  I  <l£  y  lo  ) 

The  first  and  second  examples  contained  in  Ref  15  were  run  using  FINELS. 

The  results  are  compared  with  those  obtained  by  Hinnerichs  and  Palazotto  in 
Figures  B-l  and  B-2.  The  third  example  was  not  attempted,  as  it  deals  with 
viscoelasticity  which  is  not  a  primary  concern  of  this  effort. 
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Table  VII 

Three-Bar  Linkaee  Conditions 

Param  eter 

Exam  Die  #1 

Exam  Die 

[<Tr]  o 

500  PSI 

500  PSI 

E 

30,000  PSI 

30,000  PSI 

7 

166. 667 /sec 

166.667/sec 

n 

L.O 

1.0 

F 

1  m  ax 

1250  lbs 

765  lbs 

1c 

- 

0.166667 /sec 

0 

5.0 

1.0  PSI 

dt 

3.33X10"5  sec 
_ 

3.33X10-6  sec 

Appendix  C 


Description  of  KEM Code 


The  program  FINELS  is  a  FEM  code  designed  to  predict  viscoplasticity  and 
creep  in  beam  elements  subjected  to  time-dependent  stress  or  strains.  FINELS 
can  predict  uniaxial  behavior  in  single  elements  or  can  predict  the  behavior  of 
structures  containing  several  beam  elements.  Two  versions  exist,  one  of  which 
uses  Bodner-Partom  constitutive  law  while  the  other  uses  combined  Malvern 
Overstress  theory  and  Norton's  Law. 

FINELS  calculates  elastic  and  viscoplastic  stresses  and  strains  using  the 
residual  force  method  in  which,  during  a  given  time  step  i, 

[K\{nY  =  IF]*  +  IQ}*''  (1) 

where : 

[K]  is  the  elastic  stiffness  matrix 
{u}  is  the  nodal  displacement  matrix 
[F]  is  the  matrix  of  applied  forces 
[Q]  is  the  matrix  of  residual  forces 

If  the  applied  forces  are  specified,  Equation  1  is  solved  for  {u}  by  Gaussian 
elimination.  If  the  nodal  displacements  are  specified,  ( K]  and  {u}  are  multiplied 
to  obtain  [ F] .  Each  term  Qi  of  the  residual  force  matrix  is  calculated  according 
to  the  plastic  strains  accumulated  prior  to  the  ith  time  step. 

(<?,),=  XX<r'<-  os©. 

n  —1 

(Q.)v=E£»fr'*>ne„  C2) 

n  —1 

whe  re 


N  =  To t a  1  number  of  elenients  in  system. 
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n  —Individual  element  identifier. 


0,  ^orientation  angle  of  element  n. 

The  nodal  displacements,  total  strains,  and  stresses  are  then  calculated  for 
each  member. 

{«.}•  =  ;>'! [f\ •+[<?! (3) 
=  (4) 

{*}*  =  [£]  [<  }*  _{t  -1  (5) 


Viscoplastic  and  creep  strain  rates  are  then  calculated  and  the  total  inelastic 
strain  increased  incrementally. 


{ipY= 

{i  KC}'  =  g(<r) 

{c  }•'+{«  V'C}' 


dt' 


(6) 

(?) 

(8) 


Where  the  values  of  f  and  g  are  derived  from  the  constitutive  model  or  models 
used.  Finally,  the  residual  forces  (?,■  are  recalculated  for  the  new  plastic  strain 
values . 

FINELS  has  five  parts;  a  main  program  and  four  subroutines: 

-  The  main  program  FINELS  uses  input  data  to  construct  glo¬ 
bal  and  constrained  stiffness  matrices. 

-  Subroutine  SOLVE  solves  the  matrix  equation 

\K\{u  }=[F]+[  Q\ 

for  the  nodal  displacements  {u}  by  Gaussian  elimination. 

-  Subroutine  TRN  is  the  time  integration  routine  used  when 
the  applied  forces  are  specified.  The  applied  force  i- 
assumed  to  be  a  sawtooth  wave  having  specified  maximum 


TT 


I  ill  wn1  v"o vr1  ■  vwtwvt i  nr  -w  **» '  >-%*• 


and  minimum  values  and  frequency.  Subroutine  SOLVE  is 
used  to  compute  nodal  displacements  due  to  the  applied 
forces.  Elastic  and  inelastic  stresses  and  strains  are  then 
calculated  as  previously  described  for  each  time  increment. 

-  Subroutine  DISP  is  the  time  integration  routine  used  when 
the  nodal  displacements  are  specified  as  functions  of  time. 
The  nodal  displacements  are  assumed  to  be  sawtooth  waves 
of  known  maximum  and  minimum  value  and  frequency. 
Elastic  and  inelastic  stresses  and  strains  are  then  calcu¬ 
lated  for  each  time  increment. 

-  Subroutine  READ  reads  data  describing  each  element  and 
node,  including  conditions  of  constraint. 


DATA  FILE  CREATION: 

FINELS  reads  list-directed  input  from  a  data  file  called  "input”,  which  must 
be  set  up  as  follows: 

Line:  Data  contained:  Variables: 


1 

OUTPUTl  .OUTPUT2  ,OUTPUT3 
(Output  filenames) 

OUTPUTl  =output  data  file 
0UTPUT2  =  stress-strain  plot  file 
OUTPUT3  =  time  plot  file 
(See  below) 

2 

NU  MELS,  NNODES 

NUMELS  =  No.  of  elements 
NNODES  =  No.  of  nodes 

Next  N  lines 
(N  =  NNODES ) 

X,  .  Y, 

X,-,  Y,=Xand  Y  coordinates 
of  node  i. 

Next  M 

(M  =NUMELS ) 

Ii  j  Jaj  ,Aj  ,Ej  ,It 

Ii,I2=Identific  ation  of 
first  and  second  nodes, 
A=Cross-sec  tional  area, 

E  =  Elastic  modulus, 

1-  Moment  of  inertia, 
of  element  j. 

Next  N  lines 
(N  =  NNODES ) 


Next  line 


CONSTR,,  .CONSTR,-*  .CONSTR.3  CONSTR,,  =  Co n s t r a i n t  at  node 

i  in  direc  tion  j. 

( 1  =  free  ,  O  =  constrained) 


S  or  D 


Specifies  cyclic  (S)tress 
or  ( D)isplac  e  m  e  n  t  (Strain] 


(At  this  point  the  two  models  require  different  data) 
OVERS  TRESS  MODEL  (Delete  if  Bodner  model  is  used)  : 


NextMlines  <7Y,7,n,7e  ..•?  Yield  stress,  and  viscoelastic 

(M  —  NUMELS)  and  creep  constants  for 

element  j. 

BODNER  MODEL  (De  le te  if  Ove  rstress  model  is  used)  : 


NextMlines  MINSIR 
(M=  NUMELS) 


Next  line  D0,Z,  ,Z0,Z x  ,A,m  ,n  ,r 

BOTH  MODELS  : 

Next-to-  DT,TMAX,FREQ 

last  line 

Last  line  FMAX,FMIN,UMAN,UMIN 


Minimum  stress,  below  which 
viscoplastic  strains  are  not 
calculated.  This  is  necessary 
to  avoid  dividing  by  zero 
at  very  low  stresses.  Typically 
an  arbitrary  small  number  such 
as  100  psi. 

Bodner-Partom  material  constants 


Time  increment,  total  time, 
and  frequency  of  applied 
stress  or  strain . 

Maximum  and  minimum  values 
of  applied  force  (FMAX.FMIN) 
and  displacement  (UMAX.UMIN). 


OUTPUT  FILES : 


FINELS  places  output  data  in  three  files  : 

1.  OUTPUTl .  This  file  contains  the  stress,  strain,  and  displacement  for 
each  element  during  each  time  step.  It  also  contains  the  global  and  constrained 
stiffness  matrices. 

2.  OUTPUT2  This  file  contains  the  data  necessary  to  plot  the  stress-strain 


curve  for  one  element,  arranged  in  two  columns.  The  first  column  contains 


strains  and  the  second  stresses. 

3.  OUTPUT3.  This  file  contains  the  data  requited  to  plot  stress  or  strain  vs. 
time.  If  strain  control  is  selected,  stresses  will  be  plotted.  If  stress  control  is 
selected,  strains  will  be  plotted.  The  data  is  arranged  in  four  columns,  the  first 
containing  time  and  the  rest  containing  stress  or  strain  in  elements  1,  2,  and  3 
respectively  for  each  point  in  time. 

Note  that  OUTPUT1,  OUTPLT2,  and  OUTPUT3  are  only  the  names  of  the 
internal  variables  which  identify  the  external  files.  You  can  give  these  files  what¬ 
ever  names  you  choose  in  line  1  of  "input". 

PROGRAM  FINELS : 


(See  following  pages ) 


IF  THE  BODNER- P ARTOM  MODEL  IS  USED,  USE  FINELS  AS  SHOWN.  IF 
OVERSTRESS/NORTON  MODEL  IS  USED,  SUBSTITUTE  ALTERNATE  VERSIONS 
OF  SUBROUTINES  " TRN"  AND  "DISP"  INSTEAD. 


PROGRAM  FINELS 

INTEGER  NUMELS , NNODES , DOF , II ( 10 ) , 12 ( 10 ) , Z1 , Z2 , ZMAX, N , Z , P , II , J , Q 
INTEGER  Z1MAX , Z2MAX, R , CONSTR( 40,3) 

REAL  Xl(10) ,Y1(10) ,X2(10) ,Y2(10) ,L(10) ,A(10) ,E(10) ,1(10) 

REAL  THETA(IO) ,KK( 10, 40,40) ,K( 40,40 ) ,F( 40) ,U( 40) ,K1,K2 
CHARACTER* 80  ANSWER 

OPEN( UNIT-9 , FILE- ' testin ’ , STATUS- ’ OLD ' ) 

USE  SUBROUTINE  "READ”  TO  ENTER  SYSTEM  DATA. 

CALL  READ (NUMELS, NNODES, II, 12, A, E, I, XI, Y1,X2, Y2,CONSTR,F) 
COMPUTE  LOCAL  STIFFNESS  MATRICES  FOR  GENERALIZED  BEAM  ELEMENTS 
DOF- 3 

DO  1001  N-l, NUMELS 

L(N)-SQRT( ( (X2(N)-X1(N))**2)+((Y2(N)~Y1(N) )**2)) 

THETA (N) -AS IN ( ( Y2(N)~Y1(N) )/L(N) ) 

Zl-( 3*11 ( N) )-2 
Z2-(3*I2(N) )-2 
Z1MAX— Zl+DOF— 1 
Z2MAX— Z2+DOF-1 
DO  1  Z-Zl , Z1MAX 

DO  2  P-Z , Z2MAX 

IF( P  .NE.  Z1MAX  .AND.  Z  .NE.  Z1MAX)G0  TO  3 
GO  TO  4 

IF( Z  .EQ.  Zl)GO  TO  5 
GO  TO  6 

IF( P . EQ. Z )GO  TO  7 
GO  TO  8 

K1-(A(N)*E(N) )/L(N) 

K2-12 .0*(E(N)*I(N) )/(L(N) **3 ) 

KK( N, Z , P ) -Kl* ( COS ( THETA( N) )**2)+K2*(SIN( THETA ( N ) )  **2) 
IF( P . EQ. Z+l )GO  TO  9 
GO  TO  10 

K1-(A(N) *E(N) )/L(N) 

K2-12 .0*(E(N)*I(N) )/( L( N ) *  *3 ) 
KK(N,Z,P)-(K1~K2)*SIN( THETA ( N ) ) *COS ( THETA ( N ) ) 
CONTINUE 

IF( Z . EQ. Zl+1 )GO  TO  11 
GO  TO  4 

IF( P . EQ . Z )GO  TO  13 
GO  TO  14 

K1-(A(N)*E(N) )/L( N) 

K2-12 .0*(E(N)*I(N) )/(L(N) **3 ) 

KK(N,Z,P)-(K1* ( SIN ( THETA ( N ) )**2) )  +K2* ( COS ( THETA ( N ) ) **  2 ) 


nnn 


14  CONTINUE 

4  IF( P . EQ. Z1MAX)G0  TO  15 

GO  TO  22 

15  IF( Z . EQ. Zl)GO  TO  17 

GO  TO  18 

17  KK(N,Z,P)— 6.0*( (E(N)*I(N) )/( L( N) **2 ) ) * ( SIN( THETA( N) ) ) 

18  CONTINUE 

IF( Z . EQ. Zl+1 )GO  TO  19 
GO  TO  20 

19  KK(N,Z,P)-6.0*( (E(N)*I(N) )/(L(N)**2) ) * ( COS ( THETA ( N ) ) ) 

20  CONTINUE 

IF( Z . EQ . Z1MAX)G0  TO  21 
GO  TO  22 

21  KK (N,Z,P)— 4.0*(E(N)*I(N) ) /L ( N ) 

22  CONTINUE 

KK(N,P,Z)=KK(N,Z,P) 

2  CONTINUE 

DO  23  P-Z2 , Z2MAX 
Q-P-Z2+Z1 

IF( P . LT. Z2MAX ) KK( N, Z , P )--KK( N , Z , Q) 

IF( P . EQ . Z2MAX .AND. Z . LT. Z1MAX)KK( N, Z , P )-KK( N, Z , Q) 

IF(P . EQ. Z2MAX .AND. Z . EQ. Z1MAX) KK( N, Z,P)-2.0*E(N) *1 ( N)/L( N) 

KK ( N , P , Z ) -KK ( N , Z , P ) 

23  CONTINUE 

1  CONTINUE 

DO  24  P=Z2 , Z2MAX 

DO  25  Z”P, Z2MAX 

Q-P-Z2+Z1 

R-Z-Z2+Z1 

IF( P . LT.  Z2MAX.  AND.  Z . LT.  Z2MAX) KK(N,  Z,  P )**KK(N,  R,  Q) 

IF( P. EQ. Z2MAX. AND. Z . EQ. Z2MAX) KK(N, Z, P )~KK(N, R , Q) 

IF( P .  EQ.  Z2MAX.  AND.  Z . LT.  Z2MAX)KK(N,  Z,  P)*»-KK(N/  R,  P ) 

IF( P . LT. Z2MAX. AND. Z . EQ. Z2MAX) KK(N, Z , P )“~KK( N, Z , Q) 
KK(N,P,Z)-KK(N,Z,P) 

25  CONTINUE 

24  CONTINUE 
1001  CONTINUE 

ADD  LOCAL  STIFFNESS  MATRICES  TO  GET  THE  GLOBAL  STIFFNESS  MATRIX 

P-1 
31  Z-l 

29  N— 1 

27  K(P,Z)=*K(P,Z)  +KK  ( N ,  P ,  Z ) 

N-N+l 

IF ( N . GT . NUMELS ) GO  TO  26 
GO  TO  27 

26  Z-Z+l 
ZMAX-NNODES  *DOF 

IF ( Z . GT . ZMAX ) GO  TO  28 
GO  TO  29 

28  P-P+l 

IF ( P . GT . ZMAX ) GO  TO  30 
GO  TO  31 

30  CONTINUE 

WRITE (It*) '*+*******★*****************************★★******»** ' 
WRITE ( 7, *) 'GLOBAL  STIFFNESS  MATRIX  [K] : ’ 

WRITE ( 7 , * ) ' 

DO  3100  Q-l , ZMAX 

WRITE ( 7 ,*)(K(Q,J)/J-1/ ZMAX) 


I 


3100  CONTINUE 

WRITE (7, * )'**************************************************' 
:  INPUT  APPLIED  FORCES 


3000 


Z-l 

Z-Z+l 

IF ( Z . GT . ZMAX ) GO  TO  32 

GO  TO  33 

CONTINUE 

WRITER*)  'APPLIED  FORCES  [F]  :  ' 
WRITE ( 7 , * ) ' 

DO  3000  Q-l, ZMAX 

WRITE ( 7 , * ) F(Q ) 

CONTINUE 


3200 


APPLY  DISPLACEMENT  BOUNDARY  CONDITIONS  TO  GLOBAL  STIFFNESS  MATRIX 

II-l 

J-l 

Z-(DOF*(II-l) )+J 

IF ( CONSTR ( 1 1 ,  J ) . EQ . 1 ) GO  TO  34 

GO  TO  35 

J-J+l 

IF(J.GT.DOF)GO  TO  36 

GO  TO  37 

II-II+l 

IF( II . GT . NNODES )GO  TO  38 
GO  TO  39 

F(Z)-0. 

Z1=Z 

Z2=l 

IF(Z2.EQ.Z1)G0  TO  101 
GO  TO  102 

K(Z1,Z2)-1. 

GO  TO  103 

K(Z1, Z2)-0. 

K(Z2,Z1)=0. 

Z2-Z2+1 

IF ( Z 2 . GT . ZMAX ) GO  TO  34 
GO  TO  108 
CONTINUE 

WRITE  ( 7  ,*)’**************★****★*************★*****************' 
WRITE (7, *) 'CONSTRAINED  STIFFNESS  MATRIX  [K] : ' 

WRITE (7,*)' 

DO  3200  Q= 1 , ZMAX 

WRITE ( 7 , * ) ( K( Q, J) , J-l, ZMAX) 

CONTINUE 

WRITE (7,*) ' ********************************************** • 

WRITE ( 7 , * ) 'TIME, ELEMENT, STRAIN, STRESS, DISPLACEMENT' 

CALL  THE  VISCOPLASTICITY  SUBROUTINES 


901 

2001 


READ (9,*) ANSWER 

IF ( ANSWER. EQ. 'S' )GO  TO  901 

CALL  DISP ( NUMELS , ZMAX ,X1,X2,Y1,Y2, I1,I2,L,U,E,A) 

GO  TO  2001 

CALL  TRN( NUMELS, ZMAX, F, XI, X2,Y1,Y2, I1,I2,L,K,U,E,A) 

END 

SUBROUTINE  SOLVE ( ZMAX, K, F, U) 


c 

c 

c 


DIAGONALIZE  THE  STIFFNESS  MATRIX  AND  FORCE  MATRIX 
BY  GAUSSIAN  ELIMINATION 


REAL  K( 40 , 40 ) ,F(40) ,U(40) , K1 , K2 

INTEGER  M,  N, P , Q,  ZMAX 

P-1 

304  Kl-(K(P,P)**2) 

IF(SQRT(K1) .GT. 0. 001) GO  TO  300 
GO  TO  2000 

300  M— P+1 
F(P)-F(P)/K(P,P) 

IF ( P . EQ . ZMAX ) GO  TO  305 
GO  TO  306 

305  K(P,P)-1. 0 
GO  TO  303 

306  DO  5000  N=P+1 , ZMAX 

K( P , N) =K( P , N)/K( P , P ) 

5000  CONTINUE 
K(P,P)-1.0 

302  F(M)-F(M)-( ( K  (M, P )/K( P,P))*F(P)) 

K2«K(M,P) 

DO  5001  N=P , ZMAX 

K(M,N)-K(M,N)-( (K2/K(P,P) )*K(P,N) ) 

5001  CONTINUE 
M-M+l 

IF ( M . GT . ZMAX ) GO  TO  301 
GO  TO  302 

301  P-P+l 

IF ( P . GT . ZMAX ) GO  TO  303 
GO  TO  304 

303  CONTINUE 
C 

C  SOLVE  THE  SYSTEM  OF  EQUATIONS 

C 

Q-ZMAX 
403  M-l 

401  F(M)-F(M)-((K(M,Q)/K(Q,Q))*F(Q>) 

K (M, Q ) —0 . 

M-M+l 

IF(M. GE. Q)GO  TO  400 
GO  TO  401 
400  Q-Q-l 

IF (Q. EQ. 1 )GO  TO  402 
GO  TO  403 

402  CONTINUE 
C 

C  PRINT  THE  DISPLACEMENTS 
C 

Q-l 

500  F(Q)=F(Q)/K(Q,Q) 

U(Q)=F(Q) 

Q-Q+l 

IF ( Q . GT . ZMAX ) GO  TO  2001 
GO  TO  500 

2000  PRINT*, 'CHECK  INPUTS' 

2001  RETURN 
END 

SUBROUTINE  TRN( NUMELS , ZMAX, F, XI , X2 , Y1 , Y2 , II , 12 , L, K,U , E , A) 
REAL  K(40,40),F(40),U(40),X1(20),X2(20),Y1(20),Y2(20) 
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REAL  DT, T, TMAX, STRNBP( 10 ) , EDOT( 10 ) 

REAL  FMAX, ANGLE, LI ( 10 ) , STRN1 ( 10 )  , STRAIN ( 10 ) , DISP ( 10 ) , Q 
REAL  STRESS (10) ,E(10) , A1 , A2 , L( 10 ) , KKK( 40,40), YLDSTR( 10 ) 

REAL  THETA ( 10 ) , A( 10 ) , F1P( 10 ) , F2P (10),WP(10),ZZ(10),M,NN 
INTEGER  N, 11(10) , 12 ( 10 ) , NUMELS, ZMAX, PI , P2 , P , Z 
DO  3  N-l, NUMELS 

READ  MINIMUM  STRESS 

READ( 9 , * )MINSTR( N ) 

CALCULATE  ELEMENT  ORIENTATION  ANGLES 

IF ( X2 ( N ) . GE . XI ( N ) )THETA(N)-ASIN( (Y2(N)-Y1(N) )/L(N) ) 

IF(X2(N) .LT.Xl(N) )THETA(N)-3.141593-ASIN( (Y2(N)-Y1(N) )/L(N) ) 

INITIALIZE  PARAMETERS 

STRN1 ( N) =0 . 

STRAIN(N)=0. 

STRESS(N)=0 . 

STRNBP(N)=0 . 

WP( N) =0 . 

EDOT( N) =0 . 

F1P ( N ) =0 . 

F2P ( N ) =0 . 

CONTINUE 

READ  MATERIAL  PARAMETERS,  TIME  STEP  AND  LIMIT,  FREQUENCY  (HZ), 
AND  APPLIED  MAX/MIN  FORCES  &  NODAL  DISPLACEMENTS 

READ( 9,*)DO,ZI,ZO,Z1 , AA,M, NN , R 

READ( 9 , * )DT, TMAX, FREQ 

READ ( 9 ,*) FMAX , FMIN ,  UMAX , UMIN 

DETERMINE  THE  NUMBER  OF  POINTS  TO  BE  ENTERED  AS  DATA 

QQQ-NUMELS* ( INT( TMAX/ (DT* 500 . 0 ) ) ) 

IF(QQQ.LT. 1JQQQ-1 

WRITE ( 8 , * ) STRAIN ( 1 ) , STRESS ( 1 ) 

X-1.0 
FF=0 . 

DO  4  T=0 . , TMAX, DT 

CALCULATE  APPLIED  FORCE  MATRIX 

IF( FF . GE . FMAX ) X=-l . 0 

IF(FF . LE . FMIN ) X=1 . 0 

FF-FF+ ( 2 . *FREQ* ( FMAX- FMIN ) *X*DT ) 

ANGLE =0 . 

F(l)=(FF*COS( ANGLE) ) +F1P( 1 )+FlP( 2 )+FlP( 3 ) 

F( 2 ) “ ( FF*SIN( ANGLE ) ) +F2P ( 1)+F2P(2)+F2P(3) 

DO  100  N=3 , ZMAX 
F(N)=0 . 

CONTINUE 
DO  940  P=1 , ZMAX 
DO  941  Z=1 , ZMAX 
KKK{ P , Z ) =K ( P , Z ) 

CONTINUE 

CONTINUE 


non  non  non  non  nn 


SOLVE  FOR  NODAL  DISPLACEMENTS 

CALL  SOLVE  (  ZMAX , KKK ,  F ,  U ) 

2  CONTINUE 

DO  101  N**l ,  NUMELS 

CALCULATE  ELEMENTAL  STRAINS  AND  STRESSES 

Pl-(Il(N)*3)-2 

P2-(I2(N)*3)-2 

A1-(X2(N)-X1(N)+U(P2)~U(P1) )**2 
A2-(Y2(N)-Y1(N)+U(P2+1)-U(P1+1) )**2 
LI ( N  > “SQRT ( A1+A2 ) 

STRN1 ( N)= ( Ll{ N)/L( N) ) -1 . 0 

STRESS ( N )  -  E ( N ) *  ( STRN1 ( N ) -STRNBP ( N ) ) 

IF ( STRESS (N) . GE. 0 . )XXX=1.0 
IF ( STRESS ( N ) . LT . 0 . )XXX— 1.0 

CALCULATE  VISCOPLASTIC  STRAINS 

ZZ(N)=Z1-(Z1-Z0)*EXP(-M*WP(N) ) 

IF ( ABS ( STRESS ( N ) ) . LT . YLDSTR ( N ) ) GO  TO  501 
AAA— (  (  (  ZZ ( N) /STRESS  ( N)  )  **2 . 0  )  **NN)  *  ( NN+1 . 0  )/NN 
D2P= ( DO * * 2 . 0 ) *EXP ( AAA ) 

EDOT ( N) =2 . 0*XXX* SQRT ( D2P/3 . 0 ) 

GO  TO  502 

501  EDOT ( N ) =0 . 

502  STRNBP ( N ) “STRNBP ( N )  +  ( EDOT ( N ) *DT ) 

ZREC— AA*Zl*(  ( (ZZ(N)-ZI)/Z1)**R) 

WP(N)=WP(N)+(STRESS(N)*ED0T(N)+(ZREC/(M*(Z1-ZZ(N) ) ) ) ) *DT 

CALCULATE  PLASTIC  FORCES 

F1P( N )— A( N)  *E(N) *STRNBP ( N)  *COS ( THETA (N) ) 

F2P(N)— A(N)*E(N)*STRNBP(N)*SIN(THETA(N) ) 

STRAIN ( N ) -STRN1 ( N ) 

DISP ( N) =STRAIN( N) *L(N) 

ENTER  RESULTS  INTO  OUTPUT  FILES 


Q=Q+1 

IF(Q.GE.QQQ)THEN 

Q-0 

WRITE(7,  *)'  ************************************************** 
WRITE( 7 , * )T, '1' , STRAIN ( 1 ) , STRESS ( 1) ,DISP(1) 

WRITE ( 8 ,*) STRAIN ( 1 ), STRESS ( 1 ) 

WRITE ( 10 , * ) T , STRESS ( 1 ) , STRESS ( 2 ) , STRESS ( 3 ) 

END  IF 

101  CONTINUE 

4  CONTINUE 

RETURN 
END 

SUBROUTINE  DISP ( NUMELS , ZMAX, XI , X2 , Y1 , Y2 , II , 12 , L, U, E , A ) 

REAL  XI (10) ,X2(10) ,Y1(10) , Y2 ( 10 ) , L( 10 ) ,U(10) ,E(10) ,A( 10) 
REAL  STRAIN ( 10 ) , STRNBP ( 10 ) , T, TMAX, DT, LI ( 10 ) , YLDSTR ( 10 ) , M , NN 
REAL  STRESS (10) ,K1,K2,ZZ(10) ,WP(10) , EDOT( 10 ) 

INTEGER  N, 11(10), 12(10), ZMAX, PI , P2 , NUMELS 
1-1 

DO  100  N-l, NUMELS 


read  minimum  stress 


READ( 9, *)MINSTR(N) 


INITIALIZE  PARAMETERS 


STRNBP(N)=0 . 

STRAIN(N)-0. 

STRESS(N)-0. 

WP(N)-0. 

CONTINUE 


READ  MATERIAL  PARAMETERS,  TIME  STEP  AND  LIMIT  (SEC),  FREQUENCY  (HZ), 
AND  APPLIED  MAX/MIN  FORCES  &  NODAL  DISPLACEMENTS 


READ(9,*)D0,ZI,Z0,Z1,AA,M,NN,R 
READ ( 9 , * ) DT , TMAX , FREQ 
READ ( 9 , * ) FMAX , FMI N , UMAX ,  UMIN 
WRITE ( 8 , * ) STRAIN ( 1 ) ,STRESS(1) 


DETERMINE  HOW  MANY  POINTS  TO  ENTER  INTO  OUTPUT  DATA 


QQQ-NUMELS*  ( INT(  TMAX/ ( DT* 500 . 0 )  )  ) 

IF ( QQQ . LT . 1 ) QQQ-1 

X-1.0 

UU=0. 

DO  1001  T=0 . , TMAX , DT 


CALCULATE  THE  TIME-DEPENDENT  DISPLACEMENT  MATRIX 


IF(UU.GE.UMAX)X— 1.0 

IF(UU. LE.  UMIN) X=*l .  0 

UU=UU+ ( 2 . * FREQ* ( UMAX-UMIN ) *X*DT ) 

ANGLE- 0. 

U(1)-UU* COS (ANGLE) 

U(2)=UU*SIN( ANGLE) 

DO  101  N=  3 , ZMAX 
U(N)=0. 

CONTINUE 

DO  102  N=1 , NUMELS 


CALCULATE  STRAINS  AND  STRESSES 


Pl=(Il(N)*3)-2 
P2= ( 12 (N ) *  3 ) -2 

K1-(X2(N)-X1(N)+U(P2)-U(P1) )**2 
K2-(Y2(N)-Y1(N)+U(P2+1)-U(P1+1) )**2 
LI ( N) =SQRT ( K1+K2 ) 

STRAIN ( N )  =  ( LI ( N )/L( N) )-l.  0 
STRESS(N)=E(N) * (STRAIN(N)-STRNBP(N) ) 
IF ( STRESS ( N ) . GE . 0 . )YYY=1.0 
IF(STRESS(N)  .  LT .  0  .  )YYY— 1.0 


CALCULATE  VISCOPLASTIC  STRAINS 


ZZ(N)=Z1-(Z1-Z0) *EXP( -M*WP(N) ) 

IF ( ABS ( STRESS ( N ) ) . LT . YLDSTR ( N ) ) GO  TO  501 
ZZZ—  (  (  (ZZ(N)/STRESS(N)  )  *  *2  )  *  *NN)  *  ( NN+1 .  )/NN 
D2P“(D0**2) *EXP(  ZZZ) 


ED0T(N)=2. *YYY*SQRT(D2P/3. ) 

GO  TO  502 

501  EDOT{ N) =0 . 

502  STRNBP ( N ) -STRNBP ( N )  +  ( EDOT ( N )  *DT ) 

ZREC — AA*Z1* ( ( (ZZ(N)-ZI )/Zl ) *  *R) 

XXX-ZREC/(M*(Z1-ZZ(N) ) ) 

WP(N)-WP(N)+( ( (STRESS ( N) *EDOT(N) )+XXX) *DT) 

C 

C  ENTER  RESULTS  INTO  OUTPUT  FILES 

C 

Q-Q+l 

IF ( Q . GE . QQQ ) THEN 
Q-0 

WRITE (  7 ,  * )  ’*************************************************' 
WRITE ( 7 , * ) T, N, STRAIN ( N ) , STRESS (N) ,  U(  1 ) 

IF ( N . EQ . 1 ) WRITE ( 8 , * ) STRAIN ( 1 ) , STRESS ( 1 ) 

IF(N. EQ. 1 ) WRITE ( 10 ,  * ) T/ STRESS ( 1 )  , STRESS ( 2 ) , STRESS ( 3 ) 

ENDIF 

102  CONTINUE 
1001  CONTINUE 
RETURN 
END 

SUBROUTINE  READ ( NUMELS , NNODES , II , 12 , A , E , I , XI , Y1 , X2 , Y2 , CONSTR , F ) 
REAL  A( 10 ) , E ( 10 ) , I ( 10 ) , XI ( 10 ) , X2 ( 10 ) , X( 40 ) , Y1 ( 10 ) , Y2 ( 10 ) , Y ( 40 ) 
REAL  F( 40 ) 

INTEGER  NUMELS , NNODES , II ( 10 ) , 12 ( 10 ) , CONSTR( 40 , 3 ) ,  II ,  J , Z 
CHARACTER* 20  OUTPUT1 , OUTPUT 2 , OUTPUT3 
READ( 9 , * ) OUTPUT 1 , OUTPUT2 , OUTPUT3 
OPEN(UNIT=7 , FILE=0UTPUT1 , STATUS= ' NEW' ) 

OPEN ( UNIT=8 , FILE-OUTPUT2 , STATUS- ' NEW' ) 

OPEN (UNIT- 10 , F ILE-OUTPUT 3 , STATUS- ’ NEW' ) 

C 

C  READ  NUMBER  OF  ELEMENTS  AND  NODES 

C 

READ ( 9 , * ) NUMELS , NNODES 
C 

C  READ  NODE  LOCATIONS 

C 

DO  1  11=1, NNODES 
READ (9,*)X(II),Y(II) 

1  CONTINUE 

C 

C  READ  DATA  NECESSARY  TO  DESCRIBE  EACH  ELEMENT  (NODES  AND  PROPERTIES) 
C 

DO  100  N-l, NUMELS 

READ( 9, * ) I1(N) , 12 (N) ,A(N) ,E(N) ,I(N) 

XI ( N) =X( II ( N) ) 

Y1(N)-Y(I1(N) ) 

X2 ( N) =X( 12 ( N ) ) 

Y2 ( N) =Y( 12 ( N) ) 

100  CONTINUE 
C 

C  CONDITIONS  OF  CONSTRAINT  (1  -  FREE,  0  -  CONSTRAINED) 

C 

DO  101  11=1, NNODES 

READ( 9 , * ) ( CONSTR (II,J),J=1,3) 

101  CONTINUE 
C 

C  READ  F ( X ) , F ( Y ) , AND  M(XY)  AT  EACH  NODE 


DO  102  II-l , NNODES 
Z-(3*II)-2 

READ( 9/ *)F(Z) , F ( Z+l ) , F( Z+2 ) 
102  CONTINUE 
RETURN 
END 
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★  ★ 

*  NOTE :  * 

*  -  * 

*  IF  THE  OVERSTRESS/NORTON  MODEL  IS  TO  BE  USED,  USE  THESE  * 

*  VERSIONS  OF  SUBROUTINES  " TRN"  AND  "DISP"  IN  PLACE  OF  THE  * 

*  BODNER-PARTOM  VERSIONS.  * 

*  ★ 


SUBROUTINE  TRN( NUMELS , ZMAX, F, XI , X2 , Y1 , Y2 , II , 12 , L, K, U, E, A) 
REAL  K(40,40),F(40),U(40),X1(20),X2(20),Y1(20),Y2(20) 

REAL  STRNP( 10 ) , DT , T , TMAX, YLDSTR( 10 ) , GAMMA ( 10 ) ,NEXP(10) 

REAL  FMAX , ANGLE , LI ( 10 ) , STRN1 ( 10 ) , STRAIN ( 10 ) , DISP ( 10 ) , Q 
REAL  STRESS ( 10 ) , EPDOT ( 10 ) , E( 10 ) , A1 , A2 , L( 10 ) , KKK( 40,40) 

REAL  THETA ( 10 ) , A( 10 ) , F1P ( 10 ) , F2P ( 10 ) , STRNVC( 10 ) , WP( 10 ) 

REAL  GAMC(10) ,BETA(10) 

INTEGER  N, II (10) ,12(10) , NUMELS, ZMAX, PI, P2,P,Z 
DO  3  N=l, NUMELS 

READ  YIELD  STRESS,  VISCOPLASTIC  CONSTANTS,  AND  CREEP  CONSTANTS 

READ ( 9 , * ) YLDSTR ( N ) , GAMMA ( N ) , NEXP ( N ) , G AMC ( N ) , BETA ( N ) 

CALCULATE  ELEMENT  ORIENTATION  ANGLES 

IF(X2(N) .GE.Xl(N) ) THETA ( N ) -AS IN ( (Y2(N)-Y1(N) )/L(N) ) 

IF(X2(N) .LT.Xl(N) )THETA(N)«3 . 141593~ASIN( (Y2(N)-Y1(N) )/L(N) ) 

INITIALIZE  PARAMETERS 


Q=0 

STRNl ( N ) =0 . 

STRAIN(N)=0 . 

STRNP ( N ) =0 . 

STRNVC(N)=0 . 

STRESS(N)=0 . 

F1P(N)=‘0 . 

F2P(N)-0. 

3  CONTINUE 

READ  IN  TIME  STEP  AND  LIMIT  (SEC),  FREQUENCY  (HZ),  AND  APPLIED 
MAX/MI N  FORCES  &  NODAL  DISPLACEMENTS 

READ( 9 , *)DT, TMAX, FREQ 
READ( 9 , * ) FMAX , FMIN, UMAX, UMIN 
WRITE (8,*)STRAIN(1) ,STRESS(1) 

QQQ-NUMELS* ( I NT ( TMAX/ (DT* 500 . 0) ) ) 

IF(QQQ.LT. 1)QQQ=1 

X-1.0 

FF-0. 

DO  1000  T=0 . , TMAX, DT 

CALCULATE  TIME-DEPENDENT  APPLIED  FORCE  MATRIX 
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IF(FF.GE.FMAX)X— 1.0  ! 

IF( FF . LE. FMIN )X“1 . 0  » 

FF-FF+(2.*X*FREQ*(FMAX-FMIN)*DT)  J 

ANGLE-0 .  i 

F(1)-(FF*C0S( ANGLE) ) +F1P ( 1 )+FlP( 2 )+FlP ( 3 )  j 

F{ 2 )- ( FF*SIN( ANGLE) ) +F2P ( 1 )+F2P( 2 )+F2P( 3 ) 

DO  100  N-3 , ZMAX 
F(N)-0. 

CONTINUE 
DO  940  P-1, ZMAX 
DO  941  Z-l , ZMAX 
KKK(P, Z )=K(P, Z ) 

CONTINUE 
CONTINUE 

SOLVE  FOR  NODAL  DISPLACEMENTS 

CALL  SOLVE ( ZMAX, KKK,F,U) 

DO  101  N=l, NUMELS 
AAA-0. 002 *WP(N) 

CALCULATE  STRAIN-DEPENDENT  YIELD  STRESS 

YLDSTR( N) =100000 . * ( 0 . 9+ ( 0 . 2*EXP( -AAA) ) ) 

CALCULATE  ELEMENTAL  STRAINS  AND  STRESSES 

Pl=( Il(N) *3 )-2 
P2=(I2(N)*3)-2 

A1=(X2(N)-X1(N)+U(P2)-U(P1))**2 
A2=(Y2(N)-Y1(N)+U(P2+1)-U(P1+1) )**2 
LI ( N ) =SQRT ( A1+A2 ) 

STRN1(N)=(L1(N)/L(N) )—l. 0 

STRESS ( N)=E( N) * ( STRN1 (N)-STRNP ( N)~STRNVC( N) ) 

IF( ABS ( STRESS (N) ) . LT. 10 . 0 )GO  TO  50 
GO  TO  21 
50  ZZZ=0. 

GO  TO  22 

CALCULATE  CREEP  STRAIN  RATE 

21  ZZZ=ABS( STRESS (N) )/STRESS(N) 

22  ZZQ=DT*GAMC( N )*ZZZ*( ( ABS ( STRESS ( N) ) **BETA(N) ) 

IF(ABS( STRESS (N) ) . GE . YLDSTR ( N ) ) ZZQ-0 . 

STRNVC ( N ) =ZZQ+STRNVC ( N ) 

CALCULATE  VISCOPLASTIC  STRAIN  RATE 

IF ( ABS ( STRESS ( N ) ) . GE . YLDSTR ( N ) ) GO  TO  10 
GO  TO  20 

10  ZZP-ZZZ*GAMMA(N) 

EPDOT ( N) =ZZP* ( ( ( ABS ( STRESS ( N ) )/YLDSTR(N) )-l. )**NEXP(N) ) 

WP(N) =WP(N)+( STRESS ( N) *EPDOT(N) *DT) 

STRNP ( N) = ( EPDOT ( N ) *DT) +STRNP ( N ) 

CALCULATE  NEW  PLASTIC  FORCES 

20  F1P(N)=-A(N) *E(N) * ( STRNP (N)+STRNVC(N) ) *COS(THETA(N) ) 

F2P(N) — A(N) *E(N)*( STRNP (N)+STRNVC(N) ) *SIN ( THETA( N ) ) 

STRAIN{ N) -STRN1 ( N) 


100 

941 

940 


c 

c 

c 


DISP(N)-STRAIN(N)*L(N) 


ENTER  RESULTS  INTO  OUTPUT  FILES 


101 

1000 


c 

c 

c 


c 

c 

c 


100 


c 

c 

c 

c 


c 

c 

c 


Q-Q+l 

IF ( Q . GE . QQQ ) THEN 
Q=0 

WRITE ( 7 ,*)'************************************************' 
WRITE ( 7 , * )T, N, STRAIN ( N) , STRESS( N) ,DISP(N) 

WRITE ( 7  ,  *  ) EPDOT , ZZQ , STRNP ( N )  , STRNVC ( N ) , YLDSTR ( N )  , F1P ( N ) ,F2P(N) 
WRITE ( 8 ,  * ) STRAIN ( 1 ) , STRESS ( 1 ) 

WRITE ( 10 , * ) T, STRAIN ( 1 ) , STRAIN ( 2 ) , STRAIN( 3 ) 

ENDIF 

CONTINUE 

CONTINUE 

RETURN 

END 

SUBROUTINE  DISP ( NUMELS , ZMAX, XI , X2 ) Y1 , Y2 , II , 12 , L, U, E , A ) 

REAL  XI (10) ,X2(10) ,Y1(10) ,Y2(10) ,L(10) ,U(10) ,E(10) ,A(10) 

REAL  STRAIN ( 10 ) / STRNP ( 10 ) , T, TMAX, DT, LI ( 10 ) , YLDSTR( 10 ) 

REAL  GAMMA ( 10 ) , NEXP ( 10 ) , EPDOT , STRESS ( 10 ) , K1 , K2 , STRNVC ( 10 ) 

REAL  WP(10) ,GAMC(10) , BETA(IO) 

INTEGER  N, 11(10) ,12(10) , ZMAX, PI , P2 , NUMELS 
1=1 

DO  100  N=l, NUMELS 


READ  YIELD  STRESS,  VISCOPLASTIC  CONSTANTS,  AND  CREEP  CONSTANTS 
READ ( 9 , * ) YLDSTR ( N ) , GAMMA ( N ) , NEXP ( N ) , GAMC ( N ) , BETA ( N ) 

CALCULATE  ELEMENT  ORIENTATION  ANGLES 


STRNP (N)=0. 
STRNVC ( N ) =0 . 
STRAIN ( N) =0 . 
STRESS(N)=0. 
CONTINUE 


READ  IN  TIME  STEP  AND  LIMIT  (SEC),  FREQUENCY  (HZ),  AND  APPLIED 
MAX/MIN  FORCES  &  NODAL  DISPLACEMENTS 


READ ( 9 ,*) DT , TMAX , FREQ 

READ ( 9 , * ) FMAX , FMIN , UMAX , UMIN 

WRITE ( 8 , * ) STRAIN ( 1 ) , STRESS ( 1 ) 

QQQ=NUMELS *(INT( TMAX/ (DT* 500. 0) ) ) 

IF(QQQ.LT. 1)QQQ=1 

X-1.0 

UU=0. 

DO  1001  T=0 . , TMAX,DT 


CALCULATE  TIME-DEPENDENT  DISPLACEMENT  MATRIX 


IF(UU.GE.UMAX)X=-1. 0 

IF(UU.LE. UMIN) X-1.0 

UU-UU+ ( 2 . *X*FREQ* ( UMAX-UMIN ) *DT ) 

ANGLE =0 . 

U(l)=UU*COS( ANGLE) 

U( 2)=UU*SIN( ANGLE) 

DO  101  N=3 , ZMAX 
U(N)«0. 
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1001 


CONTINUE 

DO  102  N=1 , NUMELS 

CALCULATE  ELEMENTAL  STRAINS  AND  STRESSES 

Pl=(Il(N)*3)-2 

P2-(I2(N)*3)-2 

Kl”( X2 ( N) -Xl( N) +U( P2)~U(P1) )  *  *  2 
K2”( Y2 ( N) -Y1 ( N) +U( P2+1)~U(P1+1) )  **2 
LI ( N ) -SQRT ( K1+K2 ) 

STRAIN(N)-(L1(N)/L(N) )-1.0 

STRESS( N) =E( N ) * { STRAIN ( N)~STRNP( N)~STRNVC( N) ) 
CALCULATE  CREEP  STRAIN  RATE 
ZZZ-ABS( STRESS (N) ) /STRESS (N) 

ZZQ-DT*GAMC(N) *ZZZ*( ( ABS ( STRESS( N) )**BETA(N) ) 

IF(ABS( STRESS (N) ) . GE . YLDSTR ( N ) ) ZZQ=0 . 

STRNVC ( N ) =ZZQ+STRNVC ( N ) 

CALCULATE  STRAIN-DEPENDENT  YIELD  STRESS 

AAA=0 . 002*WP(  N) 

YLDSTR (N) =100000. *(0. 9+( . 2*EXP(-AAA) ) ) 

IF ( ( ABS ( STRESS ( N ) ) ) . GE . YLDSTR ( N ) ) GO  TO  200 
GO  TO  300 

CALCULATE  VISCOPLASTIC  STRAIN  RATE 
ZZP=ZZZ*GAMMA(N) 

EPDOT=ZZP*( ( (ABS (STRESS (N) )/YLDSTR(N) )~1. )**NEXP(N) ) 
WP ( N ) =WP ( N ) + ( STRESS ( N ) *EPDOT*DT ) 

STRNP ( N ) =STRNP ( N ) + ( EPDOT*DT ) 

CONTINUE 

ENTER  RESULTS  INTO  OUTPUT  FILES 
Q-Q+l 

IF ( Q . GE . QQQ ) THEN 
Q=0 

WRITE  ( 7  ,★)'****★*****************************■******** 
WRITE ( 7 , * ) T, N, STRAIN ( N ) , STRESS ( N) , U( 1 ) 

IF( N. EQ . 1 ) WRITE ( 8 , * ) STRAIN ( 1 ) , STRESS ( 1 ) 

IF( N . EQ . 1 ) WRITE( 10 , * ) T, STRESS ( 1 ) / STRESS (2) , STRESS ( 3 ) 

END  IF 

CONTINUE 

CONTINUE 

RETURN 

END 


INPUT  FILE  "TESTIN'1  ( BODNER-PARTOM  MODEL) 


datlOO , sslOO ,  timelOO 
3,4 
0,0 
-1,1 
-1,-1 
1.4142,0 

1,2,1.0,26. 3E6 , 0 
1,3,1.0,26. 3E6 , 0 
1,4,1.0,26. 3E6 , 0 
1,1,0 
0,0,0 
0,0,0 
0,0,0 
0,0,0 
0,0,0 
0,0,0 
0,0,0 
S 

100 

100 

100 

1. E4, 600000, 915000, 1015000, 1.9E-3, 2. 57E-3, 0.7, 2. 66 
0.01,6,4. 166667E-2 
300000,-300000, .008,-. 008 
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INPUT  FILE  "TESTIN’’  (OVERSTRESS/NORTON  MODEL) 


dat413 , ss413 , time413 
3,4 
0,0 
-1,1 
-1,-1 
1.4142,0 

1, 2, 1.0,25. 05E6,0 
1,3,1.0,25. 05E6 , 0 
1,4,1. 0,25. 05E6,0 
1,1,0 
0,0,0 
0,0,0 
0,0,0 
0,0,0 
0,0,0 
0,0,0 
0,0,0 
s 

130000. . 0107.1.1.520E-26,4.022 

130000. . 0107.1.1.520E-26,4.022 

130000. . 0107.1.1. 520E-26 ,4.022 
0. 1,600,1. 666667 E~ 3 

200000. , -200000. , .008,- .008 
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and  discrepancies  noted.  The  Bodner-Partom  predictions  are  also  compared  to 
predictions  developed  using  Overstress  theory  and  Norton’s  Law  for  Secondary 
Creep  in  a  similar  FEM  code.  A  three-bar  linkage  model  is  then  developed  to 
illustrate  viscoplastic  behavior  at  a  stress  concentration.  Areas  requiring 
further  study  are  identified. 
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